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Abstract: Density matrix renormalization group (DMRG) is applied to a (1+1)- 
dimensional Xcf)'^ model. Spontaneous breakdown of discrete Z2 symmetry is studied numer- 
ically using vacuum wavefunctions. We obtain the critical coupling (A//U^)c = 59.89 ±0.01 
and the critical exponent (3 = 0.1264 ± 0.0073, which are consistent with the Monte Carlo 
and the exact results, respectively. The results are based on extrapolation to the contin- 
uum limit with lattice sizes L = 250, 500, and 1000. We show that the lattice size L = 500 
is sufficiently close to the the limit L ^ 00. 
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1. Introduction 

Hamiltonian diagonalization is a useful method for nonperturbative analysis of many-body 
quantum systems fl], • If Hamiltonian is diagonalized, the system can be analyzed non- 
perturbatively at the amplitude level using the obtained wave functions § 1^, ^, |6|, ^ . In 
addition, one can discuss associated symmetry based on operator algebra |9|. However, 
in general quantum field theories, the method does not work without reducing degrees of 
freedom because the dimension of Hamiltonian increases exponentially as the system size 
becomes large. There is a severe limitation on numerical diagonalization of Hamiltonian. 
To apply the method to quantum field theories, we need to find a way to remove unimpor- 
tant degrees of freedom and create a small number of optimum basis states. This is the 
concept of renormalization group. 

Wilson's numerical renormalization group applied to the Kondo problem was a success- 
ful consideration along this line |^0|. To analyze chain models other than Kondo Hamil- 
tonian, White proposed density matrix renormalization group (DMRG) as an extension 
of the Wilson's method [0, H. In DMRG, calculation accuracy of the target state can 



be controlled systematically using density matrices. White calculated energy spectra and 
wavefunctions of Heisenberg chains composed of more than 100 sites using a standard 
workstation. The DMRG analysis of a 100-site 5 = 1/2 chain corresponds to diagonal- 
ization of Hamiltonian with 2^'^'^ ~ 10'^'^ dimensions. DMRG has been applied to various 
one-dimensional models, such as Kondo, Hubbard, and t-J chain models, and achieved 
great success. A two-dimensional Hubbard model has also been studied with DMRG in 
both real- and momentum-space representation DMRG works well on small two- 

dimensional lattices and new techniques have been proposed for larger lattices. DMRG has 
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also been extended to finite-temperature chain models using the transfer-matrix technique 
based on the Suzuki- Trotter formula |15, 16, 17, 18, 19|. In particle physics, the massive 
Schwinger model has been studied using DMRG to confirm the well-known Coleman's pic- 
ture of 'half-asymptotic' particles at a background field = n |2^]. It is interesting to 
seek a possibility of applying the method to quantum chromodynamics (QCD) in order to 
study color confinement and spontaneous chiral symmetry breaking based on QCD vacuum 
wavefunctions. 

In fermionic lattice models, the number of particles contained in each site is limited 
because of the Pauli principle. On the other hand, in bosonic lattice models, each site can 
contain infinite number of particles in principle. It is not evident whether Hilbert space 
can be described appropriately with a finite set of basis states in bosonic models. This 
point becomes crucial when DMRG is applied to gauge theories because gauge particles 
are bosons. Before working in lattice gauge theories like Kogut-Susskind Hamiltonian pT| , 
we need to test DMRG in a simple bosonic model and recognize how many basis states 
are necessary for each site to reproduce accurate results. In this paper, we apply DMRG 
to a Xcj)'^ model with (1+1) space-time dimensions. We define a Hamiltonian model on a 
spatial lattice (only space is discretized) because DMRG is a method based on Hamiltonian 
formalism. The model has spontaneous breakdown of discrete Z2 symmetry and the exact 
values of the critical exponents are known. We are going to justify the relevance of DMRG 
truncation of Hilbert space in the bosonic model by comparing our numerical results with 



the Monte Carlo and the exact results |22, Our largest lattice size is L = 1000, which 
is about twice of the latest Monte Carlo one . It is shown that the lattice size L = 500 
is sufficiently close to the limit L — > cxd. 

This paper is organized as follows. In Sec. ^ Hamiltonian lattice formulation of the 
model is given. The canonical variables are transformed to Fock-like operators, each of 
which creates or annihilates a boson on each site. Real-space representation is maintained 
because local interactions are useful for DMRG. In Sec. ^, DMRG setup for the model is 
explained. A superblock is composed of two blocks and one site. Sec. ^ gives the definition 
of the critical values and shows numerical results. Sec. ^ is devoted to summary. 

2. A X(j)f^i Hamiltonian lattice model 

We start from the following Lagrangian of a (l+l)-dimensional Xc/)'^ model 



Lc = I dx 



(2.1) 



where space and time are continuous. Only the space is discretized to obtain a Hamiltonian 
lattice model. In this paper, the spatial derivative is modeled as a naive difference on the 
lattice. ^ Lagrangian of the lattice model is 



n=l ^ ' ^ n=l 



^Discretization errors associated with the derivative can be improved in a systematic way by introducing 
non-nearest neighbor interactions |24l ESl . 
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where a is a lattice spacing and time is continuous. Open boundary conditions are chosen. 
Hamiltonian is given by 



^ = E ii-n + + ^^^) + ^ E i<t>n - 4>n^lf , (2.3) 
n=l ^ ' ^ n=l 

where vr^ = a(j)n is the conjugate variable to All dimensionful quantities are measured 
in units of a. 

H = Ha, /lo^ = fiW, A = Aa^ (2.4) 
To quantize the system, we impose an equal-time commutation relation 

[4>m.{t),Trn{t)] = iSmn- (2.5) 



We rewrite Hamiltonian (2^) using creation and annihilation operators ali and a„. 

= (ojj + On) , vr„ = (^al - a„^ , (2.6) 



where [am, ali] = 6mn and an|0) = 0. Note that ali and are not creation and annihilation 
operators in Fock representation. The index n of the operators ali and a„ stands for the 
discretized spatial coordinate, not momentum. Real-space representation is better for our 



purpose because local interactions are useful for DMRG. (See Ref. |13, 14] for DMRG in 
momentum-space representation.) 

In this model. Fig. |l| is the only diagram that diverges 
in the continuum limit a — > 0. The divergence can be renor- 
malized by redefining the mass parameter as 



~2 ~2 

1^0 = fl -dfl , 



a 



where J/i^ is a counter term that cancels the divergence. Figure 1: The only diver 

gent diagram in the (1+1)- 
j._2 n /~2^ dimensional A^'* model. 

5"^ is given as a discrete sum because the system is quantized in a finite spatial box, 

SLil^") = wi2 , ' (2-7) 

2L ^ /~2 I A „;^2 nn 



n=i \l jl^ + 4 sin 
In the limit L ^ oo, the sum (|2.7D becomes 



7.2 



(t2 + i^)(l-t2) 



1 /tt 2 \ 
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where F is the ehiptic integral of the first kind 



F{^,k) 



de 



\/l — k? sm 



2, 



The integral Soo is used for all calculations of the counter term Sji^ even with a finite 
lattice because we are interested in the limit L ^ oo. We are going to calculate the critical 
coupling and exponent by extrapolating numerical data points to the limit L ^ oo. 

In Fock representation, the divergence can be removed easily by taking normal ordering 
of Hamiltonian. However, Hamiltonian ( |2.3D cannot be normal-ordered easily because our 
representation is given in real space and different from Fock one. For this reason, we need 
to evaluate the integral (2.8) explicitly and redefine the mass parameter to renormalize the 



divergence in a numerical manner. 



3. DMRG with one-site insertion 



We are going to apply DMRG technique to the the model given in the previous section. In 
this paper, the system is composed of two renormalized blocks (system and environment 
blocks) and one bare site, each of which is approximated with a finite number of basis 
states. Based on the DMRG technique, basis states for system and environment blocks are 
optimized to describe the whole system in a finite dimensional space. 

In fermionic models, dimension of Hamiltonian is finite on a finite lattice. On the other 
hand, in bosonic models, each site can contain any number of bosons. Namely, each site 
has infinite degrees of freedom. The dimension of a bosonic Hamiltonian is infinite even 
on a finite lattice. There is an essential difference between boson and fermion. To perform 
numerical calculations, Hamiltonian needs to be finite dimensional because a computer 
can only take care of finite dimensional matrices. For this reason, in our bosonic model, 
the number of basis states of each bare site needs to be restricted when it is inserted 
between system and environment blocks. In the White's first paper for DMRG, two sites 
are inserted between system and environment blocks in each RG step Since the 

model considered there is antiferromagnetic 5=1/2 Heisenberg chain, the total spin of 
the whole system needs to be kept constant for numerical stability. However, our model 
has no such a requirement associated with spins. When we insert bare sites between two 
blocks, a smaller number of sites is better. One-site insertion is the best choice. 

We divide Hamiltonian into two parts: 

L L-1 

^ = /^n + hn,n+l, (3.1) 
ra=l n=l 

where 

hn 

hn,n-\-l 
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We are going to apply the finite system algorithm of DMRG to the Hamiltonian. As shown 
in Fig. 0, a superblock Hamiltonian is composed of two blocks and one site: 



Hs = Hi^ + /ln-l,n + hn + h 



n,n+l 



(3.2) 



where Hi, and are effective Hamiltonian for the left and right blocks, respectively. 
hn-i,n (^n,n+i) IS an interaction between the left (right) block and the inserted n-th bare 
site. The superblock Hamiltonian is an effective Hamiltonian for the original H. 
Basis states for the superblock are 



ii,j,fc)„ = hf))iiW)i4' 



(R)\ 



(3.3) 




where 



1,...,M} and = 1,...,M} 

are finite sets of basis states for the left and right blocks, 
respectively. = 1, . . . , N} is a set of basis states 

for the n-th bare site inserted between the two blocks. 




\J 



in)\ 



(4^10), j = l,:.,N. (3.4) 




At every DMRG step, basis states ^)} and {|4^^)} 
are updated using the density matrix technique to renor- 
malize the left and right block Hamiltonian, Hi, and 
Hn. (See Appendix ^ for the creation of basis states.) 
If M and N are infinite, the superblock Hamiltonian 
Hs gives exact spectra because update of basis states 
is nothing but infinite dimensional unitary transformation. However, such a calculation is 
not possible on a finite computer. We truncate Hilbert space with finite M and N. The 
target state is expanded as 



Figure 2: DMRG sweep process is 
shown. A superblock is composed 
of two renormalized blocks (system 
and environment) and one bare site. 
A bare site is moved from end to 
end. Basis states are optimized each 
time. One DMRG sweep is com- 
posed of L basis optimization steps. 



M N M 



i=i j=i k=i 



(3.5) 



The dimension of the superblock Hamiltonian is M'^N. The relevance of truncation will 
be checked numerically by seeing convergence of energy and wavefunction with respect to 
the parameters M and A^. 

Matrix elements of the superblock Hamiltonian are 



n{i,j,k\Hs\i'j',k')n 



where 



i>n\j 



/(n)\ 
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and 



To calculate the matrix elements ( |3.6D , the following matrix elements need to be calculated 
and stored at every DMRG step. 

Other matrix elements do not need to be stored because they can be calculated with 
operator contraction. (See Appendix p.) 

4. Numerical analysis 

The critical values are defined with vacuum expGcta-tion values of the Vciriable (pn- 

m 



A _ ^ 



(4.1) 



where |\I'o) is a vacuum state and A{X) is a constant dependent on A. The vacuum expec- 
tation value V is calculated at the center n = L/2 using vacuum wavefunctions ^ijk for 
each pair of A and Jj?. The quantities /Uc(^)) /3(A), and A{X) are determined by fitting v 
to obtained numerical data. The critical coupling (A///^)c and the critical exponent (5 are 
defined in the continuum limit a ^ as follows 122 



^=(A\ +B~X, (4.2) 

/3(A) = /3 + C7A, (4.3) 

where B and C are some constants. In the continuum limit a — > 0, the coupling constant 
A = Aa^ vanishes. Then, we have 

^) =]im^, /5=lim/3(A). (4.4) 

In numerical analysis, the continuum limit cannot be reached directly. The critical values 
(A/;U^)c and /3 are determined by extrapolating A//i^(A) and /3(A) to the continuum limit 
A = 0. Also, finite-size effects on the critical values are removed by using sufficiently large 
lattices (L = 250,500, and 1000). 

Figure ^ shows N dependence of vacuum energy Eq = Eqq for each M = 8, 9, and 
10. Values of the other parameters are A = 0.6, /ig = -0.2 (/i^ = 9.744828 x 10^^), 
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and L = 1000. The energy converges as N increases with M fixed. Parameter values 
(M, N) = (10, 10) gives convergence of vacuum energy in accuracy of four digits. Figure ^ 
shows dependence of the vacuum expectation value v for each M = 8, 9, and 10. Other 
parameter values are same as Fig. ^. v converges at = 10 for each M, but convergence 
for M with N fixed is slow. We will use the values (M, N) = (10, 10) in all calculations to 
determine critical values. 

602 I 1 1 1 1 1 1 1 1 

M=8 A 
M=9 □ 
M=10 « 

601.5 - 



'0 




600.5 - 



N 

Figure 3: Vacuum energy i?o = Ei^a is plotted as a function of N for M — 8, 9, and 10. Coupling 
constant, bare mass squared, and lattice size are A = 0.6, /Iq = —0.2, and L — 1000, respectively. 



0.8 I r 

M= a A 

M= 9 □ 
0.7 -M=10 o 



V 0.4 



N 

Figure 4: Vacuum expectation value v is plotted as a function of N for AI = 8, 9, and 10. Coupling 
constant, bare mass squared, and lattice size are A = 0.6, /Iq = —0.2, and L — 1000, respectively. 



In Fig. |, vacuum expectation value v{\, /2^) is plotted as a function of A//i^ for a fixed 
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55 60 65 70 75 80 85 

Figure 5: Vacuum expectation value v is plotted as a function of X/fP with A = 0.6 fixed. The 
lattice size is L = 1000. 



A = 0.6 on the lattice L = 1000. If X/fJp is decreased from the positive side, v becomes 
smaller and finally vanishes at a point A//i^(A). The obtained data points are fitted with 
Eq. ( [4. ID to calculate /ic(A) and /5(A) for each A. The Marquardt-Levenberg nonlinear 
least-squares algorithm is used for fitting. Similar figures are drawn also for A = 1.5 and 
3.0 to determine /ic(A) and /3(A). 



(a) /j2 = -0.208 



(b| jjl = -0.2 



DMRG sweeps 



DMRG sweeps 



Figure 6: Eq (circles) and v (squares) are plotted as functions of the number of DMRG sweeps 
for A — 0.6, L — 1000, and {M,N) — (10,10). Convergence speed are compared between (a) 
/2q — —0.208 (A/2^ ~ 78.7) and (b) /i§ — —0.2 (A/i^ ^ 61.6). Circles and squares represent Eq and 
f , respectively. Data points are plotted every DMRG sweep. 



In Fig. ^, Eq and v are plotted as functions of the number of DMRG sweeps for A = 0.6 
and L = 1000. Convergence speed is compared between /2q = —0.208 (A//i^ ~ 78.7) and 
/xq = —0.2 (A//i^ ~ 61.6). The former and latter parameter values correspond to the first 
and twelfth points from the right hand side in Fig. ^, respectively. The latter is close to the 
critical value A///^(A)|^^q g. In Fig. |6|, data points are plotted every DMRG sweep. One 
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DMRG sweep is composed of L = 1000 sets of basis state optimization. Convergence tends 
to be slow as A//i^ approaches the critical value. When jj.^ = —0.208, convergence of both 
quantities is very fast. On the other hand, when /ig = —0.2, energy converges in accuracy 
of four digits even with a small number of DMRG sweeps, but convergence of v is slow, v 
barely converges in the latter case. Because of limited computer time, complete convergence 
is not available for v when is very close to the critical value f^cW- We calculate Eq and 
V with up to forty DMRG sweeps checking their convergence to determine f^ci^) on each 
lattice. 

66 I 1 1 1 1 1 r--^ „ 1 




Figure 7: The ratio A///^(A) is plotted as a function of the coupling constant A {— 0.6, 1.5,3.0) 
for three lattice sizes L = 250 (circles), 500 (triangles), and 1000 (diamonds). Each plot is fitted 
with a straight line to determine (A//i^)c in the limit A — *■ 0. The fitted lines are drawn with solid 
{L = 250), dotted {L = 500), and dot-dashed (L = 1000) lines. 



In Fig. ^, the ratio A//2^(A) is plotted as a function of A (= 0.6, 1.5, 3.0) with error bars 
for three lattice sizes L = 250, 500, and 1000. The errors come from fitting. We are going 
to take the continuum limit a — > to determine the values of {X/ ix^)c- As explained in Eq. 
( |4.2| ), the data points are fitted with a straight line and extrapolated to the limit A ^ 
for each L. The data points for A = 0.6 and 1.5 are used for fitting and extrapolation. 

Figure ^ plots the extrapolated values (A//i^)c as a function of 1/L. We observe good 
agreement between the two results for L = 500 and 1000. That is, the lattice size L = 500 
is sufficiently large and close to the limit L — > oo. For this reason, we fit the two points 
for L = 500 and 1000 with a straight line to extrapolate them to the limit L ^ oo. In 
the continuum limit a ^ and L —>■ oo, we obtain (A///^)c = 59.89 it 0.01. Our result is 



close to the Euclidean Monte Carlo result |22], which has been given with lattices up to 
L = 512. No exact result for {X/fj?)c has been known. Table [l| shows various results for 
the critical value (A//i^)c. 

Figure |9| plots the critical exponent /?(A) as a function of A (= 0.6, 1.5, 3.0) for L = 250, 
500, and 1000. As before, the data points are fitted with Eq. ( ^ ) and extrapolated to 
the limit A = for each L. The data points for A = 0.6 and 1.5 are used for fitting and 
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0.002 0.003 

1/L 



Figure 8: The critical value (A//i^)c is plotted as a function of 1/L for L = 250,500, and 1000. 
Extrapolation to the limit L ^ oo gives (A//i^)c — 59.89 ± 0.01. 



Method 



Result 



Reference 



DMRG 
Monte Carlo 

Gaussian effective potential 
Gaussian effective potential 
Connected Green function 
Coupled cluster expansion 
Non-Gaussian variational 
Discretized light cone 
Discretized light cone 



59.89 ±0.01 



61.56 



+0.48 
-0.24 



61.266 
61.632 
58.704 
22.8 < {\/n'^)c < 51.6 
41.28 
43.896, 33.000 
42.948, 46.26 



This work 

m 
m 



m 

m 

m 



Table 1: Various results for the critical coupling constant (A//i^)c are listed. Our result is consistent 
with Monte Carlo within the errors. 



extrapolation. Figure |l^ plots the extrapolated values (3(0) as a function of 1/L. Also in 
this case, we observe that the results for L = 500 and 1000 are very close. Extrapolation 
with these two points to the limit L ^ oo gives the critical exponent /3 = 0.1264 it 0.0073. 
No Monte Carlo result for the critical exponent /3 is available in the literature at this 
point. According to the universality class consideration, the critical exponents of the 
(l+l)-dimensional X(j)^ model are same as the two-dimensional Ising model, which has 



been exactly solved by Onsager [23|. Our result is consistent with the exact value /3 = 
1/8 = 0.125. As seen from Fig. ^ and 10, the lattice size L = 500 is sufficiently close to 
the limit L — > oo. 



5. Summary 

We have studied spontaneous breakdown of Z2 symmetry of a (l+l)-dimensional cp'^ model 
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A 

Figure 9: /3(A) is plotted as a function of the coupling constant A (= 0.6, 1.5,3.0) for three lattice 
sizes L = 250 (circles), 500 (triangles), and 1000 (diamonds). Each plot is fitted with a straight line 
to determine /3(0). Fitted lines are drawn with solid (L=250), dotted {L = 500), and dot-dashed 
(L = 1000) lines. 

0.15 
0.145 

0.14 
0.135 

0.13 

f](0) 0.125 
0.12 
0.115 
0.11 
0.105 
0.1 

0.001 0.002 0.003 0.004 0.005 

1/L 

Figure 10: (3{0) is plotted as a function of 1/L for L = 250,500, and 1000. Extrapolation to the 
limit Z/ — > oo gives (3 = 0.1264 ± 0.0073, which is consistent with the exact result /3 = 1/8 = 0.125. 

using the density matrix renormalization group technique. We have determined the critical 
coupling constant (A//x^)c and the critical exponent /? of the model by extrapolating the 
numerical results for finite lattices to the continuum limit a ^ and L — oo. DMRG 
truncation works well also in the bosonic model. The lattice with L = 500 can give results 
sufficiently close to the limit L ^ oo. To improve calculations near ficWi we need to find 
a way to include effects of large quantum fluctuations. 
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A. Basis states and density matrices 

This is a brief explanation of basis state creation witli the density matrix technique, which 
is the key ingredient of DMRG. We split a spatial lattice into left and right blocks, each of 
which are expressed with basis sets = 1, 2, . . . , iVL} and = 1,2,..., N^}, 

respectively. The target state to be solved is expanded as 

I*)=EE*^^-I^^'^^)l^'^''^)' ^L<A^R. (A.l) 

i=l j=l 

When A'^L and N-^ are sufficiently small, we can obtain wavefunction of a target state (say 
ground state) by diagonalizing Hamiltonian numerically. If the basis states and 
are transformed, the wavefunction is also transformed. We want to make the 
absolute values of the wavefunction components very small as many as possible, because 
basis states giving very small wavefunction component are not important for the target 
state 1^) and can be thrown away. We transform basis states based on the following 
singular-value decomposition. 

^ij = Y,UikDkVkj, (A.2) 

k=l 



where 

In the new basis 
the target state becomes 



i=i j=i 

A^L ATr 

1=1 j=i 



ATl 



m=J2Dk\nt%r)- (A.3) 

k=l 

The result does not change largely even if a basis state with small is removed from the 
calculation. On the other hand, a basis state with large Dk is important and cannot be 
neglected. We can use D^'s to choose good basis states and control calculation accuracy. 
In actual numerical works, we diagonalize the following density matrices to obtain U, V, 
and D in stead of performing singular-value decomposition directly. 

Nn ATl 

= T.'^h^^'j = T.^:k\Dk\'Ui'k, 
j=i k=i 

A^L ATl 

= j:n'^ij' = J2Vkj\DkfVkr. 

i=l k=l 
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See Ref. [. 



19 1 for the details of DMRG and related topics. 



B. Matrix elements 



This appendix gives formulas for matrix elements of powers of the field operators 7r„ and 
4>n- For the definition of the basis states Ij*-""^), see Eq. (|3.4| ). 



V(j-l)(j-2)<5,_i,,.+i 



+(2j - l)Sjj, - V(j'-l)(i'-2)(^,+i,,v_i 



1 

1 



V(j-l)0--2)5,_iy+i 



+(2j - l)S,y + V(/-l)(i'-2)5,+i,v_i 



V(j-l)(i-2)(j-3)(j-4)<5,_i,,.+3 



+4jV(i - 1)(J - 2)5. 



-6V(j-l)(j-2)5,_i,,v+i 
+3(2j2 - 2i + 1)5,- ,v 



+4jV(j'-l)(/-2)5i+2j' 
-6v'(j'-l)(/-2)5,+i,,'_i 

+ \/(/ - W - 2)(j' - 3)(j' - 4)5,+3,i' 
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